function pca_analysis()
% PCA_ANALYSIS 主成分分析专项实验
% 
% 深入研究基于SVD的主成分分析方法和应用

fprintf('=== 主成分分析专项实验 ===\n\n');

%% 实验1：PCA基础理解
fprintf('实验1：PCA基础理解\n');
fprintf('目标：理解PCA的数学原理和几何意义\n\n');

% 生成2D相关数据
n_samples = 200;
rng(42);  % 固定随机种子

% 原始数据（在主轴方向有不同方差）
X_original = randn(n_samples, 2);
X_original(:, 1) = 3 * X_original(:, 1);  % x方向方差更大
X_original(:, 2) = 1 * X_original(:, 2);  % y方向方差较小

% 旋转数据
theta = pi/4;  % 45度旋转
R = [cos(theta), -sin(theta); sin(theta), cos(theta)];
X1 = (R * X_original')';

fprintf('数据集信息:\n');
fprintf('样本数量: %d\n', n_samples);
fprintf('特征维度: %d\n', size(X1, 2));
fprintf('数据均值: [%.4f, %.4f]\n', mean(X1));
fprintf('数据标准差: [%.4f, %.4f]\n', std(X1));

% 数据中心化
X1_centered = X1 - mean(X1);
fprintf('中心化后均值: [%.4f, %.4f]\n', mean(X1_centered));

% 协方差矩阵
C1 = cov(X1_centered);
fprintf('协方差矩阵:\n');
disp(C1);

% 使用SVD进行PCA
[U1, S1, V1] = svd(X1_centered, 0);
principal_components = V1;
explained_variance = diag(S1).^2 / (n_samples - 1);

fprintf('主成分分析结果:\n');
fprintf('第一主成分方向: [%.4f, %.4f]\n', principal_components(:, 1));
fprintf('第二主成分方向: [%.4f, %.4f]\n', principal_components(:, 2));
fprintf('解释方差: [%.4f, %.4f]\n', explained_variance);
fprintf('方差解释比例: [%.2f%%, %.2f%%]\n', ...
        100 * explained_variance / sum(explained_variance));

% 验证：与特征值分解比较
[eigenvectors, eigenvalues] = eig(C1);
[eigenvalues_sorted, idx] = sort(diag(eigenvalues), 'descend');
eigenvectors_sorted = eigenvectors(:, idx);

fprintf('\n特征值分解验证:\n');
fprintf('特征值: [%.4f, %.4f]\n', eigenvalues_sorted);
fprintf('特征向量1: [%.4f, %.4f]\n', eigenvectors_sorted(:, 1));
fprintf('特征向量2: [%.4f, %.4f]\n', eigenvectors_sorted(:, 2));

% 可视化
figure('Name', 'PCA基础理解', 'Position', [100, 100, 1200, 400]);

subplot(1, 3, 1);
scatter(X1(:, 1), X1(:, 2), 'b.', 'MarkerSize', 8);
axis equal; grid on; title('原始数据');
xlabel('X1'); ylabel('X2');

subplot(1, 3, 2);
scatter(X1_centered(:, 1), X1_centered(:, 2), 'b.', 'MarkerSize', 8);
hold on;
% 绘制主成分方向
pc1_line = 3*sqrt(explained_variance(1)) * [-principal_components(:, 1)'; principal_components(:, 1)'];
pc2_line = 3*sqrt(explained_variance(2)) * [-principal_components(:, 2)'; principal_components(:, 2)'];
plot(pc1_line(:, 1), pc1_line(:, 2), 'r-', 'LineWidth', 3);
plot(pc2_line(:, 1), pc2_line(:, 2), 'g-', 'LineWidth', 2);
axis equal; grid on; title('中心化数据和主成分');
legend('数据点', 'PC1', 'PC2');

% 主成分空间投影
X1_pca = X1_centered * principal_components;
subplot(1, 3, 3);
scatter(X1_pca(:, 1), X1_pca(:, 2), 'r.', 'MarkerSize', 8);
axis equal; grid on; title('主成分空间');
xlabel('PC1'); ylabel('PC2');

%% 实验2：高维数据PCA
fprintf('\n实验2：高维数据PCA\n');
fprintf('目标：在高维数据上应用PCA降维\n\n');

% 生成高维数据
n_samples2 = 150;
n_features = 10;
true_dim = 3;  % 真实内在维度

% 生成低维潜在数据
Z = randn(n_samples2, true_dim);
% 生成随机投影矩阵
W = randn(n_features, true_dim);
% 高维观测数据
X2 = Z * W' + 0.2 * randn(n_samples2, n_features);

fprintf('高维数据信息:\n');
fprintf('样本数量: %d\n', n_samples2);
fprintf('观测维度: %d\n', n_features);
fprintf('真实内在维度: %d\n', true_dim);

% 数据标准化
X2_std = (X2 - mean(X2)) ./ std(X2);

% PCA分析
[U2, S2, V2] = svd(X2_std, 0);
explained_variance2 = diag(S2).^2 / (n_samples2 - 1);
explained_variance_ratio = explained_variance2 / sum(explained_variance2);
cumulative_variance = cumsum(explained_variance_ratio);

fprintf('主成分分析结果:\n');
fprintf('成分    特征值    方差解释比    累积解释比\n');
fprintf('----    ------    ----------    ----------\n');
for i = 1:min(8, length(explained_variance2))
    fprintf('%4d    %6.3f    %10.3f    %10.3f\n', ...
            i, explained_variance2(i), explained_variance_ratio(i), cumulative_variance(i));
end

% 选择主成分数量
variance_threshold = 0.95;
n_components = find(cumulative_variance >= variance_threshold, 1);
fprintf('保留%.0f%%方差需要 %d 个主成分\n', variance_threshold*100, n_components);

% 降维和重构
X2_reduced = U2(:, 1:n_components) * S2(1:n_components, 1:n_components);
X2_reconstructed = X2_reduced * V2(:, 1:n_components)';

% 重构误差
reconstruction_error = norm(X2_std - X2_reconstructed, 'fro') / norm(X2_std, 'fro');
fprintf('重构相对误差: %.4f\n', reconstruction_error);

%% 实验3：PCA在数据可视化中的应用
fprintf('\n实验3：PCA在数据可视化中的应用\n');
fprintf('目标：使用PCA将高维数据可视化\n\n');

% 生成多类别高维数据
n_classes = 4;
n_samples_per_class = 50;
n_features3 = 8;

X3 = [];
labels = [];
class_centers = randn(n_classes, n_features3) * 3;

for c = 1:n_classes
    % 每个类别的数据围绕不同中心分布
    class_data = randn(n_samples_per_class, n_features3) + class_centers(c, :);
    X3 = [X3; class_data];
    labels = [labels; c * ones(n_samples_per_class, 1)];
end

fprintf('多类别数据信息:\n');
fprintf('类别数: %d\n', n_classes);
fprintf('每类样本数: %d\n', n_samples_per_class);
fprintf('总样本数: %d\n', size(X3, 1));
fprintf('特征维度: %d\n', n_features3);

% 标准化数据
X3_std = (X3 - mean(X3)) ./ std(X3);

% PCA降维到2D
[U3, S3, V3] = svd(X3_std, 0);
X3_pca_2d = U3(:, 1:2) * S3(1:2, 1:2);

% 计算2D投影的方差解释比例
explained_var_2d = sum(diag(S3(1:2, 1:2)).^2) / sum(diag(S3).^2);
fprintf('2D投影保留方差比例: %.2f%%\n', explained_var_2d * 100);

% 可视化
figure('Name', 'PCA数据可视化');
colors = ['r', 'g', 'b', 'm'];
markers = ['o', 's', '^', 'd'];

for c = 1:n_classes
    class_mask = labels == c;
    scatter(X3_pca_2d(class_mask, 1), X3_pca_2d(class_mask, 2), ...
            50, colors(c), markers(c), 'filled');
    hold on;
end

grid on; axis equal;
xlabel('第一主成分'); ylabel('第二主成分');
title(sprintf('PCA 2D可视化 (保留%.1f%%方差)', explained_var_2d * 100));
legend('类别1', '类别2', '类别3', '类别4');

%% 实验4：PCA用于特征选择
fprintf('\n实验4：PCA用于特征选择\n');
fprintf('目标：使用PCA进行特征选择和降维\n\n');

% 生成具有不同重要性的特征
n_samples4 = 100;
n_important_features = 3;
n_noise_features = 7;
n_total_features = n_important_features + n_noise_features;

% 重要特征（有结构的数据）
X_important = randn(n_samples4, n_important_features);
X_important(:, 2) = X_important(:, 1) + 0.5 * randn(n_samples4, 1);  % 相关特征
X_important(:, 3) = 2 * X_important(:, 1) + 0.3 * randn(n_samples4, 1);  % 强相关特征

% 噪声特征（随机数据）
X_noise = 0.1 * randn(n_samples4, n_noise_features);

% 合并数据
X4 = [X_important, X_noise];

fprintf('特征选择数据信息:\n');
fprintf('总特征数: %d\n', n_total_features);
fprintf('重要特征数: %d\n', n_important_features);
fprintf('噪声特征数: %d\n', n_noise_features);

% 标准化
X4_std = (X4 - mean(X4)) ./ std(X4);

% PCA分析
[U4, S4, V4] = svd(X4_std, 0);
explained_variance4 = diag(S4).^2 / (n_samples4 - 1);
explained_variance_ratio4 = explained_variance4 / sum(explained_variance4);

% 分析主成分载荷
loadings = V4 * S4 / sqrt(n_samples4 - 1);

fprintf('主成分载荷分析:\n');
fprintf('特征    PC1载荷    PC2载荷    PC3载荷    重要性得分\n');
fprintf('----    -------    -------    -------    ----------\n');

for i = 1:n_total_features
    importance_score = sum(loadings(i, 1:3).^2);  % 前3个主成分的载荷平方和
    fprintf('%4d    %7.3f    %7.3f    %7.3f    %10.3f\n', ...
            i, loadings(i, 1), loadings(i, 2), loadings(i, 3), importance_score);
end

% 特征重要性排序
[~, feature_importance_order] = sort(sum(loadings(:, 1:3).^2, 2), 'descend');
fprintf('特征重要性排序: [%s]\n', num2str(feature_importance_order'));

%% 实验5：增量PCA
fprintf('\n实验5：增量PCA\n');
fprintf('目标：实现和测试增量PCA算法\n\n');

% 生成大数据集（分批处理）
total_samples = 1000;
batch_size = 100;
n_features5 = 5;
n_batches = total_samples / batch_size;

fprintf('增量PCA设置:\n');
fprintf('总样本数: %d\n', total_samples);
fprintf('批次大小: %d\n', batch_size);
fprintf('批次数量: %d\n', n_batches);

% 初始化增量PCA参数
mean_incremental = zeros(1, n_features5);
cov_incremental = zeros(n_features5, n_features5);
n_seen = 0;

% 生成和处理数据批次
fprintf('批次处理结果:\n');
fprintf('批次    累积样本数    前3个特征值\n');
fprintf('----    ----------    -----------\n');

for batch = 1:n_batches
    % 生成当前批次数据
    X_batch = randn(batch_size, n_features5);
    % 添加一些结构
    X_batch(:, 2) = X_batch(:, 1) + 0.5 * randn(batch_size, 1);
    X_batch(:, 3) = 0.3 * X_batch(:, 1) + 0.7 * X_batch(:, 2) + 0.2 * randn(batch_size, 1);
    
    % 更新均值
    batch_mean = mean(X_batch);
    new_mean = (n_seen * mean_incremental + batch_size * batch_mean) / (n_seen + batch_size);
    
    % 更新协方差矩阵（简化版本）
    X_batch_centered = X_batch - batch_mean;
    batch_cov = (X_batch_centered' * X_batch_centered) / (batch_size - 1);
    
    if n_seen == 0
        cov_incremental = batch_cov;
    else
        % 简化的协方差更新
        cov_incremental = (n_seen * cov_incremental + batch_size * batch_cov) / (n_seen + batch_size);
    end
    
    mean_incremental = new_mean;
    n_seen = n_seen + batch_size;
    
    % 计算当前的主成分
    [~, eigenvals] = eig(cov_incremental);
    eigenvals = sort(diag(eigenvals), 'descend');
    
    fprintf('%4d    %10d    [%.3f, %.3f, %.3f]\n', ...
            batch, n_seen, eigenvals(1:3));
end

% 与批量PCA比较
fprintf('\n增量PCA vs 批量PCA比较:\n');
% 生成完整数据集进行比较
X_full = randn(total_samples, n_features5);
X_full(:, 2) = X_full(:, 1) + 0.5 * randn(total_samples, 1);
X_full(:, 3) = 0.3 * X_full(:, 1) + 0.7 * X_full(:, 2) + 0.2 * randn(total_samples, 1);

% 批量PCA
X_full_centered = X_full - mean(X_full);
[~, S_full, ~] = svd(X_full_centered, 0);
eigenvals_batch = diag(S_full).^2 / (total_samples - 1);

fprintf('方法        前3个特征值\n');
fprintf('----        -----------\n');
fprintf('增量PCA     [%.3f, %.3f, %.3f]\n', eigenvals(1:3));
fprintf('批量PCA     [%.3f, %.3f, %.3f]\n', eigenvals_batch(1:3));
fprintf('相对误差    [%.2e, %.2e, %.2e]\n', ...
        abs(eigenvals(1:3) - eigenvals_batch(1:3)) ./ eigenvals_batch(1:3));

%% 实验6：核PCA概念演示
fprintf('\n实验6：核PCA概念演示\n');
fprintf('目标：演示非线性数据的PCA局限性\n\n');

% 生成非线性数据（螺旋形）
n_samples6 = 200;
t = linspace(0, 4*pi, n_samples6);
spiral_data = [t .* cos(t); t .* sin(t)]' + 0.5 * randn(n_samples6, 2);

fprintf('非线性数据信息:\n');
fprintf('数据类型: 螺旋形\n');
fprintf('样本数量: %d\n', n_samples6);
fprintf('原始维度: 2\n');

% 标准PCA
spiral_centered = spiral_data - mean(spiral_data);
[U_spiral, S_spiral, V_spiral] = svd(spiral_centered, 0);
explained_var_spiral = diag(S_spiral).^2 / (n_samples6 - 1);
explained_ratio_spiral = explained_var_spiral / sum(explained_var_spiral);

fprintf('线性PCA结果:\n');
fprintf('第一主成分解释方差比例: %.2f%%\n', explained_ratio_spiral(1) * 100);
fprintf('第二主成分解释方差比例: %.2f%%\n', explained_ratio_spiral(2) * 100);

% 简单的多项式特征扩展（模拟核方法）
X_poly = [spiral_data, spiral_data.^2, spiral_data(:,1).*spiral_data(:,2)];
fprintf('多项式特征扩展后维度: %d\n', size(X_poly, 2));

% 扩展特征上的PCA
X_poly_centered = X_poly - mean(X_poly);
[U_poly, S_poly, V_poly] = svd(X_poly_centered, 0);
explained_var_poly = diag(S_poly).^2 / (n_samples6 - 1);
explained_ratio_poly = explained_var_poly / sum(explained_var_poly);

fprintf('多项式核PCA结果:\n');
fprintf('前2个主成分解释方差比例: %.2f%%\n', sum(explained_ratio_poly(1:2)) * 100);

% 可视化比较
figure('Name', '线性PCA vs 核PCA概念', 'Position', [100, 100, 1200, 400]);

subplot(1, 3, 1);
scatter(spiral_data(:, 1), spiral_data(:, 2), 'b.', 'MarkerSize', 8);
axis equal; grid on; title('原始螺旋数据');

subplot(1, 3, 2);
spiral_pca = spiral_centered * V_spiral;
scatter(spiral_pca(:, 1), spiral_pca(:, 2), 'r.', 'MarkerSize', 8);
axis equal; grid on; title('线性PCA投影');

subplot(1, 3, 3);
poly_pca = X_poly_centered * V_poly;
scatter(poly_pca(:, 1), poly_pca(:, 2), 'g.', 'MarkerSize', 8);
axis equal; grid on; title('多项式特征PCA投影');

%% 实验7：PCA在图像处理中的应用
fprintf('\n实验7：PCA在图像处理中的应用\n');
fprintf('目标：使用PCA进行图像特征提取\n\n');

% 生成简单的"人脸"图像数据集
img_size = 32;
n_images = 50;

% 创建基础人脸模板
[X_img, Y_img] = meshgrid(1:img_size, 1:img_size);
base_face = exp(-((X_img - img_size/2).^2 + (Y_img - img_size/2).^2) / (img_size/3)^2);

% 生成变化的人脸图像
face_images = zeros(n_images, img_size^2);
for i = 1:n_images
    % 添加随机变形
    shift_x = randn() * 2;
    shift_y = randn() * 2;
    scale = 1 + randn() * 0.1;
    
    face = exp(-(((X_img - img_size/2 - shift_x).^2 + (Y_img - img_size/2 - shift_y).^2) / (img_size/3)^2)) * scale;
    face = face + 0.1 * randn(img_size, img_size);  % 添加噪声
    
    face_images(i, :) = face(:);
end

fprintf('人脸图像数据集:\n');
fprintf('图像数量: %d\n', n_images);
fprintf('图像尺寸: %d×%d\n', img_size, img_size);
fprintf('特征维度: %d\n', img_size^2);

% 图像PCA（特征脸）
face_mean = mean(face_images);
face_centered = face_images - face_mean;

[U_face, S_face, V_face] = svd(face_centered, 0);
eigenfaces = V_face;  % 特征脸
explained_var_face = diag(S_face).^2 / (n_images - 1);
explained_ratio_face = explained_var_face / sum(explained_var_face);

fprintf('特征脸分析:\n');
fprintf('前5个主成分解释方差比例:\n');
for i = 1:5
    fprintf('  PC%d: %.2f%%\n', i, explained_ratio_face(i) * 100);
end

cumulative_face = cumsum(explained_ratio_face);
n_components_90 = find(cumulative_face >= 0.9, 1);
fprintf('保留90%%方差需要 %d 个特征脸\n', n_components_90);

% 图像重构演示
test_image_idx = 1;
original_image = reshape(face_images(test_image_idx, :), img_size, img_size);

reconstruction_ranks = [1, 3, 5, 10, 20];
fprintf('\n图像重构质量:\n');
fprintf('主成分数    重构误差    PSNR(dB)\n');
fprintf('--------    --------    --------\n');

for r = reconstruction_ranks
    % 投影到主成分空间
    coefficients = face_centered(test_image_idx, :) * eigenfaces(:, 1:r);
    % 重构
    reconstructed = face_mean + coefficients * eigenfaces(:, 1:r)';
    reconstructed_image = reshape(reconstructed, img_size, img_size);
    
    % 计算重构误差
    mse = mean((original_image(:) - reconstructed_image(:)).^2);
    psnr = 10 * log10(1 / mse);
    reconstruction_error = norm(original_image(:) - reconstructed_image(:)) / norm(original_image(:));
    
    fprintf('%8d    %8.4f    %8.2f\n', r, reconstruction_error, psnr);
end

fprintf('\n=== PCA分析实验总结 ===\n');
fprintf('1. PCA通过SVD实现，提供最优线性降维\n');
fprintf('2. 主成分数量选择需要平衡信息保留和维度减少\n');
fprintf('3. PCA在数据可视化中非常有效\n');
fprintf('4. 主成分载荷可以用于特征选择\n');
fprintf('5. 增量PCA适合大规模数据处理\n');
fprintf('6. 线性PCA对非线性结构有局限性\n');
fprintf('7. PCA在图像处理中有广泛应用\n\n');

end